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2. DEVELOPMENT OF AN EXPLICIT TIME MARCHING PROCEDURE 
FOR LAMINAR AND TURBULENT FLOW 

- SUMMARY VIEWGRAPHS 
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AN EXPLICIT FINITE-VOLUME TIME -MARCHING PROCEDURE 
FOR TURBULENT FLOW CALCULATIONS 

Start: Denton explicit time-marching method. 
Allure - easy to understand method. 

Continuity 

dp + dpu + dpv + dpw = 0 
at dx ay az 

6p = [- apu - dpv - dpw ] <5t 
dx dy dz 

Momentum 

6{pu) = [steady eqn momentum error] 6t 
<S(pv) = [ ] St 

6(pw) = [ ] St 
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Start: Denton explicit time-marching method. 

Questions : 

1. Is smoothing necessary for convergence of explicit method ? 


2. Why, at low Mach numbers, is the CFL criterion used to get the 
time step for the momentum equations ? 

( 6t = <Sx/ [velocity + speed of sound] ) 

3. Why not extend the method to laminar and turbulent flow ? 

What are the problems involved ? 

4. Why does he use an interpolated pressure in the momentum 
equation for transonic flow ? 

Can we show why and when it is stable ? 

5. How can the method be extended to separated flow ? 
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ANSWER > 

Development of Explicit method for calculation of 
Xnviscid, Laminar or Turbulent Flow 
Mach number = 0 to >2.5, including shocks 
Economical - grid points 
With or without separation 
Tested on 2-d duct flows 
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1. Is smoothing necessary for convergence of explicit method ? 
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New" control volume 
3 unknowns • 

3 (we 11 -posed) equations 


YES 


NO 


( "New" control volume, traditionally used for boundary layers ) 
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2. Why, at low Mach numbers, is the CFL criterion used 
to get the time step for the momentum equations ? 

( St = 6x/ [velocity + speed of sound] ) 

CONSERVATIVE FORM OF MOMENTUM EQUATION 

d(pu) + V’pii u = - dp + . . . 

at ax 

U [ dp 4- V’pu ] + pdU + Pii'Vu = -ap + . . . 

at at ax 

continuity 

included, therefore <5t . = 6t 

cont mom 

Stability analysis, continuity and momentum > 

CFL condition 6t s 6x/(u+c) 

CONVECTIVE FORM OF EQUATION 

pau + pu*Vu = -ap + . . . 
at ax 

Stability analysis, momentum equation > <St = <5x/u 
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3. Why not extend the method to laminar and turbulent flow ? 
What are the problems involved ? 

RESOLUTION OF TRANSVERSE PRESSURE GRADIENT 



Flat plate turbulent boundary layer ap/ay n 0 

P = pRT 

6p = RT 6p - pRU 6U 

6p dependent on continuity and momentum errors - 
stability is highly grid and <5t dependent 
difficult without smoothing 

Borrow idea from pressure correction methods - 
6p depends oftly on continuity error. 

6p = RT 6p 

= [ -apu - apv - dpw ] RT 6t 
ax ay az 

Stable without smoothing. Multi-volume approach 
needed for highly nonuniform 6y grid spacing. 
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4. Why does Denton use an interpolated pressure in 
the momentum equation for transonic flow ? 

Can we show why and when it is stable ? 


Want apu + apv + apw = 0 

dx ay az 

1-D stability analysis. p = p + 6p 

u r u + 6u 
6u - - C a(6p)/ax 
6p 'z RT ap 

Interpolated pressure 


P i =(pRT) i or P i =. 5[(pRT) i+1 +(pRT) i _ 1 ] or . . . 


> 


A i 6p i 


+ A 


i+l 6p i+l + A i-l 6p i-l 
= continuity error. 


+ . . . 

each control volume 


Explicit method approximation 

[l/6t] <Sp ^ = continuity error for control volume 


Stability requires positive and dominant. 
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5. How can the method be extended to separated flow ? 
(a) UPWIND DIFFERENCING 


> u o • o 

li-1 I i i+1 


u > 0 u du/dx = u (u i - u i _ 1 )/6x 

U < 0 u du/dx = (-u)(u i - u i+1 )/6x 

> positive coefficient for u. 

i 

(b) UPWINDED CONTROL VOLUMES 

- control volumes depend on local u 
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4. BACKFLOW - EXTENSIONS TO THE COMPUTATIONAL PROCEDURE 
4a. Discretization of Convection Terms 


The momentum and energy equations are discretized over control 
volumes fixed relative to the grid points. Central differencing is 
used except in regions where there are large cross flows or 
backflow. In these regions a side upwind or reverse upwind 
differencing is used for stability. The details are as follows. 

Control volume for momentum or energy for point i+l,j. 

X x j + 1 

N 


W x x E j 


S 

X X j-1 

i i + 1 


Bulk flow 
direction 
> 


Convection of property <$ where 
d> = u for x-momentum 

= v for y-momentum 

* = h for energy (enthalpy) equation. 

Convection term integrated over control volume 

/pvi’Vd) dVol = /V*pu<fc dVol - (fc^/V-pu dVol 

= £4>PiA‘dA - 4> M ^PlA'dA 

= (pu*A) N <t> N + (PU’A) s <t> s + (PU’A) E <t> E + (PU'A) w <t> w 

- 4> M t(PirA) N + (Pil'A)g + (pu*A) e + (pu*A) w 3 

We wish to express this in terms of the <$>’s at the grid points, 
i.e we want the equation in the form 


/pu'Vtfc dVol 


C E*i+l,j + C W*i,j + C EE*i+2,j 


The coefficients C are determined from the mass fluxes through the 
sides (pu*A) and the discretization choice for <J> N , d» s , 4>g, $ W3 and 


For stability we wish the center point coefficient, CL,, to be 

Ei 

positive and greater than the sum of the other positive 
coefficients. 



When [(PU'A) N +(PU-A) s +(PU-A) E +(PirA) w ] > 0 


take <J> M = ^ . 

to give a negative contribution to C^. 
When t(py-A) N +(PU-A) s +(PU*A) E +(Pll-A) w ] < 0 


take + M =*i + l,o 

to give a positive contribution to Cg. 

<t> E and (j> w 


When ( pu * A ) E > 0 
*E = *i+l,j 

This centered evaluation of <J>g (second order accurate) gives a 
positive contribution to Cg. 

When (pu-A) e < 0 

*E = ^i+2,0 

This upwind evaluation of <{>g (first order accurate) gives a 
negative contribution to Cgg. 

This also determines <J>^ since <J>g for one control volume is dfy for 
the next control volume. 
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and <J>g 


The geometrical centered evaluation for is 

*N = 0.5 [<1-F)(* i+1< + F(* i+1 ,j + l+* i ,j + l> 1 

where F is the fraction of the distance of the North face between 
the grid points. 


j+1 


-N - 


F = a/b 


For accuracy this centered evaluation should be used whenever 
possible. 

For stability when 


take 


(pu*A) n > 0 and (Pii*A)g > 0 


F ^ (pm’A) e /(P1A-A) n 


When the inequality is chosen, which for equal grid spacing will 
occur when (pu*A) n > 2(pia*A) e » 

(pu*A) N <t> N = + <fr i( j) C(Pil - A) N - (PU’A) E ] 

+ 0 - 5( '"i + l,j + l + *i,j + l ,<PU ' &) E 
When (pu*A) n > 0 and < Pii ’ A) E 1 0 

(the primary flow is locally backwards or zero relative to the 
bulk flow direction) 

take 

By symmetry for (pu’^)g > 0 and (pu'A) e > 0 


take 


F 2 1 - (PU'A^/tPU'A}. 




For (pu*A) s > 0 and (pia'A) e 0 


and 


take 


*S = 


When (pu*A)^ is negative, (pii'A)g is positive for the next, the 
j + 1, control volume, so that is determined from <t>g for the j+1 
control volume. Similarly when (pu*A)g is negative, <fcg is 
determined from for the j-1 control volume. 

Comparison with earlier scheme. 


In these terms, Nicholson (Section 3, Report 5, JM/86-6) 
considered only positive values for (pu'A)g, i.e. no reverse flow. 

The formulae he used for <J>^, and were the same as given 

here. However the upwinding he took for the cross flows was 
different. In particular when (pu-£)^ was positive, <J)^ was 

evaluated using 

F i 1 <Pu-A) E(at 1 _ 1 ) /(pu-4) n + 1 

Taking F i, ( PU*A) E /( PU*A) N gives lower and hence more conservative 
values of F when (pu‘A) N > 2(Pii*A) E , i.e., where for uniform grid 
spacing, the geometric F may not be used. 



high low 

cross flow cross flow 


F used in calculations. 
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Table 1 SIDE UPWINDING ALGORITHM 
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4b. Improved Pressure Interpolation for SBLI 


The same Mach number dependent pressure interpolation formula 
for the calculation of the density used earlier is also used here. 
However to converge calculations with strong shock boundary layer 
interactions, i.e. with shock induced separations, the Mach number 
used in the formula needed to be changed from the local Mach 
number to the local free stream Mach number. Since stability is 
compromised if the Mach number is underestimated but not if it is 
overestimated, for simplicity, the value used was the largest Mach 
number on the relevant pair of i-surfaces. 





M for these points 
max 


P i + 1 = ^ (p i + a o (p i + i- p i ) * I 1 (p i + r p i-i J 

a 0 (M) + aj(M) + a 2 <M) = 1 
M = Mach number = maximum Mach number at 


+ I 2 < p i + l- p i-2> ! 

planes i and i+1. 
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4c. Evaluation of Turbulent Viscosity for the Present Test Cases 


The turbulence model used in the calculations is a Prandtl 
mixing length formulation 

= pL 2 "du/dy” 

where the mixing length L is the smaller of 

0. 41y (with a Van Driest correction) 
or 0.086. 

Sometimes it is difficult to determine an appropriate boundary 
layer thickness, 6. 

For the present test cases 

(a) UTRC, boundary layer separation geometrically triggered and 

(b) Sajben diffuser, separation induced by shock, 

the boundary layer thickness used to calculate L was changed to 
obtain a reasonable separation when compared with the 
measurements. The details of what was used to determine the 
effective boundary layer thickness follow. 

Determining the edge of the boun dary layer. 

The location of the edge of the boundary layer is determined 
by the total pressure gradient Idp^/dyl. In particular a search 

starts from outside the boundary layer (in the middle of the duct) 
and proceeds towards the wall to locate where 

(p - p) *DPFACT 
Q • 

(local duct width) 

This is the edge of the boundary layer for the mixing length 
calculation. The larger the DPFACT the smaller the boundary layer 
thickness . 



Case DPFACT 

(a) UTRC, sep. b.l. 275 

(b) Sajben, p /p = 0.722 2.5 

6 O 
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"Time" lag for boundary layer thickness. 


For case (b), Sajben p /p = 0.722, shrinking 6 by increasing 

© o 

DPFACT was insufficient to correctly obtain the separation induced 
by the shock. Qualitatively since turbulence is convected with the 
flow there needs to be time for the turbulence to change - it does 
not change suddenly. This was qualitatively introduced into the 
calculation by lagging the boundary layer thickness used for the 
calculation of L by 5 grid points. (The lag is between 0.5 and 
0.75 throat heights through the separation region. ) In particular 
after i = 40 (x/h=1.7, upstream of the shock at x/h~2.4 but well 
downstream of the throat at x/h=0) the mixing length in the outer 
part of the boundary layer was obtained using 

L(i) = 0. 086( i-5 ) . 

The time lag was used only for case (b). 
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5. BACKFLOW - TEST CASES 


The extensions to the computational procedure described in 
section 4 were necessary for modelling two extreme cases of 
separated flow, the UTRC separated and reattached flat plate 
turbulent boundary layer (NASA Contract NAS3-22770, reference 1) 
and the MDRL transonic diffuser flow with a strong shock (MDRL 
Report No. 81-07, reference 2). These cases exhibit large boundary 
layer blockage (displacement thickness/local duct height), large 
backflow velocities, relative to the free stream velocity, and 
high rms/mean turbulence levels in the backflow region. The 
maximum boundary layer blockages were 58 percent (fifty eight!) in 
the UTRC low speed (U rg ^ = 27 m/s) flow and 27 percent in the MDRL 

diffuser with a shock Mach number of 1.353. The maximum backflow 
velocities were 37 percent and 25 percent, respectively, of the 
local maximum free-stream velocity. The ratios of rms/mean axial 
velocities at the locations of maximum reverse flow velocity were 
35 percent in the UTRC flow and 66 percent in the MDRL diffuser. 

If the backflows in the two cases were varying sinusoidally, these 
values would correspond to maximum backflow velocities of -5.4 + 
2.7 m/s amd -71.7 +67.3 m/s, respectively. 

UTRC Separated and Reattached Turbulent Boundary Layer 


The geometry and streamlines for flow through the UTRC test 
section are shown in Fig. 1; and laser doppler velocity 
measurements are shown in Fig. 2. The corresponding calculated 
velocity vectors together with the locus of points for which U=0 
are seen in Fig. 3. The size of the reverse flow region is well 
modelled, and the maximum calculated reverse flow velocity of -4.1 
m/s agrees well with the measured maximum value of -5.4 m/s. This 
good agreement for the reverse flow leads to reasonable agreement 
between the measured and calculated values of skin friction 
coefficient in the separation zone, as shown in Fig. 4. The 
calculated locations of separation and reattachment are seen to be 
close to the measured locations. The good agreement in the 
separated flow region was obtained with the Prandtl mixing length 
model by reducing the turbulent viscosity in the boundary layer as 
discussed in section 4 (i.e. by using DPFACT = 2.0). This then 
gave a corresponding decrease in the calculated skin friction 
upstream and downstream of the separation zone, as seen in Fig. 4. 
We conclude that the present explicit computational procedure can 
be used for flows with extensive and strong backflow but that a 
more sophisticated turbulence model is required. 
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MDRL Diffuser - Strong Shock Case 


With a back pressure, P exit /P Q inlet * of 0.722, the MDRL diffuser 

G had a shock Mach number of 1.353. Shock induced separation 
occurred in the turbulent boundary layer on the curved top wall. 
This contrasts with the case of 0.805 pressure ratio which gave a 
shock Mach number of 1.235 and no separation. In this section, 
results of calculations for these two flows will be compared, with 
particular attention being given to the backflow in the strong 
shock case. 

The calculated shock locations are clearly seen for the two 
cases in the contours of static pressure in Fig. 5. The strong 
shock is located further downstream and shows evidence of a lambda 
foot at the curved top wall. The computed shocks are both quite 
sharp as a result of the use of the M&M pressure interpolation 
formula (see section 3 of this report, reference 3). 

For the strong shock case, the computed and measured static 
pressure distributions on the top wall are compared in Fig. 6. The 
computed shock is just downstream of the measured 1 coat ion and is 
therefore somewhat stronger with a shock Mach number of 1.39. 
Upstream of the shock the static pressures are indistinguishable; 
but downstream the calculated static pressures are consistently 
higher than those measured, perhaps partly because of 
three-dimensionality in the measured flow. 

The Mach number contours in Fig. 7 show the flow accelerating 
up to the shock and decelerating downstream. The top wall boundary 
layer thickens appreciably more through the strong shock. This is 
seen also in the velocity vectors of Fig. 8, which show the 
separation bubble downstream of the strong shock. Fig. 9 shows 
this calculated backflow in more detail, and for comparison the 
magnitude and possible variations of the measured backflow are 
also shown. The maximum calculated backflow velocity of -87.7 m/s 
agrees quite well with the maximum measured value of -71.7 m/s. 

Figures 8 and 9 demonstrate quite graphically the significant 
blockage caused by the separation bubble; and this is also seen in 
Fig. 10, which shows contours of total pressure. 

We conclude that calculations of diffuser flows with strong 
shocks and shock induced separation can be performed with the 
present explicit method. As discussed in section 4, this 
calculation required a time lag of the turbulent viscosity to give 
a reduced viscosity in the separation bubble. In fact, this simple 
modification to the turbulence model produced a dramatic upstream 
movement of the shock and the calculation rapidly converged on a 
shock location close to that measured. Again this suggests the 
need for a more sophisticated turbulence model. But the present 
study of strong backflows has clearly demonstrated that they can 
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be modelled with an explicit method based on the finite volume 
approach . 
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Fig. 3 Calculated velocity vectors and the locus of points for which U 








LOCAL STATIC/INLET TOTAL PRESSURE 



Fig. 6 Computed and measured static pressure distributions 
on the curved top wall of the MDRL diffuser. 

Strong shock case; P exit /P 0inlet = °* 722 * 
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Fig. 7 SAJBEN TRANSONIC DIFFUSER Mach Number, Contour Interval = 0.05 



.805 




Calculated velocity vectors upstream of and within the separation bubble 





Fig. 10 SAJBEN TRANSONIC DIFFUSER Total Pressure , P T /P n , Contour Interval = 0.025 


6. MACH NUMBER DEPENDENT INTERPOLATION FORMULA 
FOR DENSITY-UPDATE TIME-MARCHING METHODS 


A 1-d stability analysis of density-pressure relations used in 
the computation of transonic flow was performed in Report No. 
JM/85-11 (see section 3 of this report, reference 3). Here we give 
a parallel development of a density interpolation equation for 
effective pressure for use in density-update methods. The formulae 
considered are tested using the density-update scheme outlined in 
Table 1. 

Downwind Effective Pressure 


In section 2.8 of reference 3, we considered an inconsistency 
in the pressure-density relation such that the pressure used in 
the momentum equation is offset by one grid point from the density 
used in the continuity equation, i.e. 

P i = P i+1 RT - (40) 

In a density update method this may be viewed as an effective 
pressure evaluated downwind of its point of use in the momentum 
equations. This pressure-density relation was found to be stable 
for all Mach numbers, but it results in poor shock capturing as 
the calculated shock is spread over numerous grid points. Fig. 1 
shows the calculated and theoretical pressure distributions for a 
1-d calculation with a nominal shock number of 1.45. 


Mach Number Dependent Interpolation Formula for Effective Pressure 


In section 2.6 of reference 3, we saw that when the Mach 
number is high, the density update method is stable with the ideal 
gas equation of state satisfied at each grid point, i.e.. 


p . = p . RT . 
l li 


(73) 


Since this is the correct pressure-density relation for ideal 
gases it should be used where feasible. In this section, we will 
start with a generalized density interpolation equation for 
effective pressure 


P 1 = p i + l + a o <p i- p i + l> + I 1 (p r p i + 2 ) + (p i- p i + 3> 


(74) 


with 


p i = 


e. RT. 
P i i 


(75) 
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and 


1 


(76) 


a + a 1 + a 0 = 
o 1 2 

for second order accuracy. 

We seek Mach number limitations to a Q , a^ and a^ using the 
stability criterion that 

the center point coefficient must be greater than 
the sum of the other positive coefficients, 

Coef center > Sum Coef + ' <46) 

Substituting 


f£i = < VlL + f2 ) ' 5p i + a ' a o >6p i + l 


- ®l sp i+2 


" i2 6p i+3 


into Eq. 25 of reference 3 and rearranging in terms of the 
coefficients of each 6p^, aQ, a^, and yields (neglecting 

variations of A, u and c with i) 



Now let us consider a simple second order scheme with a 0 = 0 and 
a^ = 1 - a^, and find limiting values of a Q . From Eq. 74, it is 
obvious that we should consider only values in the range 
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(79) 


0 1 a Q i i. 


In this range, the coefficient of <s P ^ + 3 is positive or zero, and 
for the coefficient of 6p^ + ^ (the center point) to be greater than 
the coeffcient of <Sp^ + g, we require 


a <2+2 M«(M+1). (80) 

° 5 5 

The coefficient of 6p^ is positive when 

2a - M*(M+1) > 0 or M(M+1) < 2a /S. (81) 

o o 

In this region, from Eq. 46, we require 


2 + M*(M+1) - 3a > 1 

o o 


a + 2a - MS(M+1) 
_Q o 

2 


or 


a < 1 + 4 M*(M+1) . 

° 3 9 


(82) 

(83) 


This criterion is more restrictive that Eq. 80 and the 
corresponding stability limit is shown as a function of Mach 
number in Fig. 2 for the conservative case of S = 1.0. 

A set of equations for a Q , a^ and a which satisfy the 

stability criteria (Eqs. 80 and 83) and give second order accurate 
interpolation (Eq. 76) has been selected; that is 

a =4 M(M+1 ) 

° 9 

a. = 1 - a (84) 

1 o 



0 . 


This Mach number dependent formulation for a Q and is shown in 

Fig. 3. These equations are referred to as the M&M Mach number 
dependent interpolation formula for density-update time-marching 
methods . 
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Computational Tests of the M&M Density Interpolation Formula 


In this section, results of shock capturing with the M&M 
formula (Eq. 84) in the density update method (Table 1) are 
presented for Denton’s 1-d nozzle. 

Calcula tion Details 

Number of axial grid points 
At inlet i 

For air k 

^exit^t, inlet 
Results 

The variations of static pressure, Mach number, and total 
pressure for all three back pressures are shown in Fig. 4. All 
three shocks are captured over four steps. The upstream side of 
the shock is sharply defined with only minor deviations from the 
theoretical 1-d solution. On the downstream side, there is a small 
overshoot and undershoot in static pressure and Mach number over 
two steps; the total pressure distributions show no overshoots or 
undershoots amd show a sharp decrease over two steps. 

Concluding Remarks 

It is hoped that the M&M density interpolation formula will be 
useful to those organizations like NASA Lewis who are using 
density-update time-marching codes. It is also hoped that the 
stability analysis performed under this NASA Grant will be 
enlightening to users and developers of time-marching codes. 


= 46, 6x = 1 

= 1, M = 0.80 

=1.4, R = 287. J/kgK 

= 0.85, 0.80, 0.75 
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Table 1. 


Outline of Effective Pressure Method with Different Time Steps 
Density Update Scheme. 

UNKNOWNS ( 2 - DIMENSIONS) 

^,u,v,y?u) ,^>v) , (^e o ),h o ,P,T 


CONTINUITY 

J^r + V»ym) - 0 

%f> ■ (- ^ yu)« dA) ^t/Vol 

/ Y + ^ 

p - ^?RT 

MOMENTUM (IN VISCID) 
^_(^u) + (/;u) u = -V p 


Cl) 


( 2 ) 


USING EQ. 2 - (u or v) TIMES EQ. 1 

y?&u « (-^ (uy?u)«d4 + pi^dA) 

+ (j) u)»dA))*§ t^/Vol 

v » (- ^ (v^>u) *dA + pjl»dA) 

+ <f") * dA) )* S t m /V°l 

u * u + Su 

v » v + $ v 

y>u) - f>u 
(f>v) -j> v 


ENERGY 


h * constant 
o 

T - (h. - u 2 +v 2 )/C_ 
o -j” P 
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X 

Fi g. 1 Comparison of calculated and theoretical 1-D static 

pressure distributions, PW = P/P . _ 
r t, inlet 

theoretical; 

calculated using a downwind effective pressure, Eq. 4 

Grid spacing, Sx = 1; p exit / p t>inlet = °* 8 0* 
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STABILITY LIMITS, A0+A1=1 


M (M + 1) 


0.00 0.20 0.40 0.80 0.80 1.00 1.20 1.40 1.80 1.80 2.00 2.20 2.40 2.60 2.80 3.00 

MACH 


Fig. 2 M & M density interpolation formula for use with 

Effective Pressure, Density Update, Time Marching Methods. 

Stability limits for a Q when a 2 - 0 and a Q + - 1. 




MflCH NUMBER DEPENDENT PI’S WITH R0+R1=1 



0.00 0.20 0.00 0.S0 0.80 1. 00 1.20 1.00 1.80 1.80 2.00 2.20 2.00 2.80 2.00 3.00 

MRCH 


Fig. 3 M & M Mach number dependent values for the density interpolation 

coefficients, a^ and a^, in the Effective Pressure Method. 
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DENTON ID EXAMPLE AO=M* (M+l ) *4/9 Al=l-AO 


pw = 


p/p 


t, inlet 



o 


0.0 10.0 20.0 30.0 > 40.0 50.0 


X 


Fig. 4a Calculated 1-D solution for Denton's nozzle 
using M & M density interpolation formula 
with the Nicholson/Moore Effective Pressure, 
Density Update Method. 


Calculations for three exit static pressures at x = 46., 
P exit /P t, inlet ' °' 85 ' °- 80 " and °' 75 - 
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DENTON ID EX 



0.860 0.880 0.900 0.920 0.940 0.960 0.980 1.000 



